
# ------ Libraries ------ #
library(texreg)
library(xergm)
set.seed(1912)

setwd("YOUR WORKING DIRECTORY")



# --- Notification:

note <- function(message=message) {system(sprintf("/usr/local/bin/terminal-notifier -title 'RStudio' -message '%s' -activate org.rstudio.RStudio",
                                                  message), ignore.stdout=TRUE, ignore.stderr=TRUE, wait=FALSE)}




# --- Data:

mid.nw_full <- readRDS("data/pauls_cranmer/prepped/system/mid.nw_full.rds")
directcont.nw_full <- readRDS("data/pauls_cranmer/prepped/system/directcont.nw_full.rds")
capratio_full <- readRDS("data/pauls_cranmer/prepped/system/capratio_full.rds")
tradedep_full <- readRDS("data/pauls_cranmer/prepped/system/tradedep_full.rds")
igosec_full <- readRDS("data/pauls_cranmer/prepped/system/igosec_full.rds")
igoecon_full <- readRDS("data/pauls_cranmer/prepped/system/igoecon_full.rds")



# ------ Models ------ #

# --- Model 1 (Baseline):

model.1_baseline <- btergm(
  mid.nw_full ~ 
    edges + 
    altkstar(2, fixed = TRUE) + 
    cycle(4) + 
    gwesp(0, fixed = TRUE) + 
    nodemix("democ", -3) + 
    edgecov(directcont.nw_full) + 
    edgecov(capratio_full) + 
    edgecov(tradedep_full) + 
    edgecov(igosec_full) + 
    edgecov(igoecon_full) + 
    memory(type = "autoregression", lag = 1), 
  R = 2000,
  parallel = "multicore", ncpus = 3);note("Done.")

screenreg(model.1_baseline)

saveRDS(model.1_baseline, file = "inference/system/model.1_baseline.rds")











# --- Model 2 (Baseline + Weak Comms):
  
dir <- "data/comm_detection/detected_comms/weak/"
file_list <- list.files("data/comm_detection/detected_comms/weak/", pattern = "weak_ties_")


for(i in 1:length(file_list)){

  file.i <- file_list[[i]]
  
  print(i)
  
  weak_comms <- readRDS(paste(dir, file.i, sep = ""))[1970:1989] 
  
  model.2_weak <- btergm(
    mid.nw_full ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(directcont.nw_full) + 
      edgecov(capratio_full) + 
      edgecov(tradedep_full) + 
      edgecov(igosec_full) + 
      edgecov(igoecon_full) + 
      edgecov(weak_comms) +
      memory(type = "autoregression", lag = 1), 
    R = 2000,
    parallel = "multicore", ncpus = 3);note("Done.")
  
  name <- paste("model.2_weak", file.i, sep = "_")
  
  saveRDS(model.2_weak, file = paste("inference/system/weak/", name, sep = ""))

  
}




# --- Average across different settings and save down 

m1 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_50d_15x.rds")))
m2 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_50d_25x.rds")))
m3 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_100d_15x.rds")))
m4 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_100d_25x.rds")))
m5 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_200d_15x.rds")))
m6 <- as.data.frame(summary(readRDS("inference/system/weak/model.2_weak_weak_ties_200d_25x.rds")))

model2_avg <- (m1 + m2 + m3 + m4 + m5 + m6)/6

model2_avg_tr <- createTexreg(coef.names = rownames(model2_avg), 
                        coef = model2_avg$Estimate, 
                        ci.low = model2_avg$`2.5%`,
                        ci.up = model2_avg$`97.5%`)

screenreg(model2_avg_tr)
#saveRDS(model2_avg_tr, file = "models/model2_avg_tr.rds")












# --- Model 3 (Baseline + Strong Comms):

dir <- "data/comm_detection/detected_comms/strong/"
file_list <- list.files("data/comm_detection/detected_comms/strong/", pattern = "strong_ties_masternodes_")

for(i in 1:length(file_list)){
  
  file.i <- file_list[[i]]
  
  print(i)
  
  strong_comms <- readRDS(paste(dir, file.i, sep = ""))[1971:1990] 
  
  model.3_strong <- btergm(
    mid.nw_full ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(directcont.nw_full) + 
      edgecov(capratio_full) + 
      edgecov(tradedep_full) + 
      edgecov(igosec_full) + 
      edgecov(igoecon_full) + 
      edgecov(strong_comms) +
      memory(type = "autoregression", lag = 1), 
    R = 2000,
    parallel = "multicore", ncpus = 3);note("Done.")
  
  name <- paste("model.3_strong", file.i, sep = "_")
  
  saveRDS(model.3_strong, file = paste("inference/system/strong/", name, sep = ""))
  
}




# --- Average across different settings and save down 

m1 <- as.data.frame(summary(readRDS("inference/system/strong/model.3_strong_strong_ties_masternodes_prop0.2.rds")))
m2 <- as.data.frame(summary(readRDS("inference/system/strong/model.3_strong_strong_ties_masternodes_prop0.25.rds")))
m3 <- as.data.frame(summary(readRDS("inference/system/strong/model.3_strong_strong_ties_masternodes_prop0.3.rds")))

model3_avg <- (m1 + m2 + m3)/3 

model3_avg_tr <- createTexreg(coef.names = rownames(model3_avg), 
                              coef = model3_avg$Estimate, 
                              ci.low = model3_avg$`2.5%`,
                              ci.up = model3_avg$`97.5%`)

screenreg(model3_avg_tr)
#saveRDS(model3_avg_tr, file = "models/model3_avg_tr.rds")











# --- Model 4 (Baseline + Strong Comms -- no contig):

dir <- "data/comm_detection/detected_comms/strong/"
file_list <- list.files("data/comm_detection/detected_comms/strong/", pattern = "strong_ties_masternodes_")

for(i in 1:length(file_list)){
  
  file.i <- file_list[[i]]
  
  print(i)
  
  strong_comms <- readRDS(paste(dir, file.i, sep = ""))[1971:1990] 
  
  
model.4_strong_nocontig <- btergm(
  mid.nw_full ~ 
    edges + 
    altkstar(2, fixed = TRUE) + 
    cycle(4) + 
    gwesp(0, fixed = TRUE) + 
    nodemix("democ", -3) + 
    edgecov(capratio_full) + 
    edgecov(tradedep_full) + 
    edgecov(igosec_full) + 
    edgecov(igoecon_full) +
    edgecov(strong_comms) +
    memory(type = "autoregression", lag = 1), 
  R = 2000,
  parallel = "multicore", ncpus = 3);note("Done.")


name <- paste("model.4_strong_nocontig", file.i, sep = "_")

saveRDS(model.4_strong_nocontig, file = paste("inference/system/strong/", name, sep = ""))

}




# --- Average across different settings and save down 

m1 <- as.data.frame(summary(readRDS("inference/system/strong/model.4_strong_nocontig_strong_ties_masternodes_prop0.2.rds")))
m2 <- as.data.frame(summary(readRDS("inference/system/strong/model.4_strong_nocontig_strong_ties_masternodes_prop0.25.rds")))
m3 <- as.data.frame(summary(readRDS("inference/system/strong/model.4_strong_nocontig_strong_ties_masternodes_prop0.3.rds")))

model4_avg <- (m1 + m2 + m3)/3 

model4_avg_tr <- createTexreg(coef.names = rownames(model4_avg), 
                              coef = model4_avg$Estimate, 
                              ci.low = model4_avg$`2.5%`,
                              ci.up = model4_avg$`97.5%`)

screenreg(model4_avg_tr)
#saveRDS(model4_avg_tr, file = "models/model4_avg_tr.rds")












# --- Model 5 (Network effects + node effects for weak comms):

dir <- "data/pauls_cranmer/prepped/subsets/weak_nodes/outcome_nets/"
file_list <- list.files(dir, pattern = "mid_sub_w_")

directcont_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/directcont_sub_w.rds")
capratio_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/capratio_sub_w.rds")
tradedep_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/tradedep_sub_w.rds")
igosec_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igosec_sub_w.rds")
igoecon_sub_w <- readRDS("data/pauls_cranmer/prepped/subsets/weak_nodes/predictors/igoecon_sub_w.rds")

for(i in 1:length(file_list)){

  file.i <- file_list[[i]]
  
  print(i)
  
  mid_sub_w <- readRDS(paste(dir, file.i, sep = ""))
  

weak_sub_full <- btergm(
  mid_sub_w ~ 
    edges + 
    altkstar(2, fixed = TRUE) + 
    cycle(4) + 
    gwesp(0, fixed = TRUE) + 
    nodemix("democ", -3) + 
    edgecov(directcont_sub_w) +
    edgecov(capratio_sub_w) + 
    edgecov(tradedep_sub_w) + 
    edgecov(igosec_sub_w) + 
    edgecov(igoecon_sub_w) +
    nodematch("weak_comm") +
    nodefactor("weak_bridge") +
    memory(type = "autoregression", lag = 1), 
  R = 2000,
  parallel = "multicore", ncpus = 3);note("Done.")


name <- paste("model.5_mid_sub_weak", file.i, sep = "_")

saveRDS(weak_sub_full, file = paste("inference/community_subset/weak/", name, sep = ""))

}





# --- Average across different settings and save down 

m1 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_50d_15x.rds")))
m2 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_50d_25x.rds")))
m3 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_100d_15x.rds")))
m4 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_100d_25x.rds")))
m5 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_200d_15x.rds")))
m6 <- as.data.frame(summary(readRDS("inference/community_subset/weak/model.5_mid_sub_weak_mid_sub_w_200d_25x.rds")))

model5_avg <- (m1 + m2 + m3 + m4 + m5 + m6)/6

model5_avg_tr <- createTexreg(coef.names = rownames(model5_avg), 
                              coef = model5_avg$Estimate, 
                              ci.low = model5_avg$`2.5%`,
                              ci.up = model5_avg$`97.5%`)

screenreg(model5_avg_tr)
#saveRDS(model5_avg_tr, file = "models/model5_avg_tr.rds")












# --- Model 6 (Network effects + node effects for strong comms):

dir <- "data/pauls_cranmer/prepped/subsets/strong_nodes/outcome_nets/"
file_list <- list.files(dir, pattern = "mid_sub_s_weaknodes_")

directcont_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/directcont_sub_s_weaknodes.rds")
capratio_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/capratio_sub_s_weaknodes.rds")
tradedep_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/tradedep_sub_s_weaknodes.rds")
igosec_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igosec_sub_s_weaknodes.rds")
igoecon_sub_s <- readRDS("data/pauls_cranmer/prepped/subsets/strong_nodes/predictors/nodeset_weak/igoecon_sub_s_weaknodes.rds")


for(i in 1:length(file_list)){

  file.i <- file_list[[i]]
  
  print(i)
  
  mid_sub_s <- readRDS(paste(dir, file.i, sep = ""))
  
  
  strong_sub_full <- btergm(
    mid_sub_s ~ 
      edges + 
      altkstar(2, fixed = TRUE) + 
      cycle(4) + 
      gwesp(0, fixed = TRUE) + 
      nodemix("democ", -3) + 
      edgecov(directcont_sub_s) +
      edgecov(capratio_sub_s) + 
      edgecov(tradedep_sub_s) + 
      edgecov(igosec_sub_s) + 
      edgecov(igoecon_sub_s) +
      nodematch("strong_comm") +
      nodefactor("strong_bridge") +
      memory(type = "autoregression", lag = 1), 
    R = 2000,
    parallel = "multicore", ncpus = 3);note("Done.")
  
  
  name <- paste("model.6_mid_sub_strong", file.i, sep = "_")
  
  saveRDS(strong_sub_full, file = paste("inference/community_subset/strong/", name, sep = ""))
  
}




# --- Average across different settings and save down 

m1 <- as.data.frame(summary(readRDS("inference/community_subset/strong/model.6_mid_sub_strong_mid_sub_s_weaknodes_prop0.2.rds")))
m2 <- as.data.frame(summary(readRDS("inference/community_subset/strong/model.6_mid_sub_strong_mid_sub_s_weaknodes_prop0.25.rds")))
m3 <- as.data.frame(summary(readRDS("inference/community_subset/strong/model.6_mid_sub_strong_mid_sub_s_weaknodes_prop0.3.rds")))

model6_avg <- (m1 + m2 + m3)/3

model6_avg_tr <- createTexreg(coef.names = rownames(model6_avg), 
                              coef = model6_avg$Estimate, 
                              ci.low = model6_avg$`2.5%`,
                              ci.up = model6_avg$`97.5%`)

screenreg(model6_avg_tr)
#saveRDS(model6_avg_tr, file = "models/model6_avg_tr.rds")


